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LOW THRUST SPACE VEHICLE TRAJECTORY 
OPTIMIZATION USING REGULARIZED VARIABLES 

INTRODUCTION 


The space vehicle trajectory optimization problem in general is treated from the 
point of view of the Mayer-Bolza problem which is well known in the calculus of 
variations [1,2]. Usually the solution must be obtained numerically, requiring some 
means of iteration procedures involving multiple integrations of sets of differential 
equations. One of the primary considerations in evaluating a numerical optimization 
procedure is the computation time for a predetermined accuracy with which the terminal 
conditions have to be satisfied. It is known that, when space vehicle trajectories cross 
regions with strongly varying gravitational force fields, the necessary' numerical accuracy 
often requires extreme computation times. The computation time depends on ( 1 ) the 
numerical integration procedure, (2) the iterative numerical optimization method, and 
(3) the mathematical formulation of the problem used. 

In this study, standard integration [3] and standard iteration [4,5] techniques are 
used. The main purpose of the investigation was to develop an improved formulation of 
the set of differential equations describing the space vehicle motion and the optimality 
conditions; in other words, to find a set of equations with high stability of numerical 
integration and less sensitivity with respect to errors in the required guesses for the 
unknown boundary conditions. These characteristics are necessary for the iteration 
procedure to demonstrate good convergence characteristics in isolating the optimal 
solution. 

In celestial mechanics, regularizing transformations are used to remove the 
singularity which occurs due to the r" 2 term in the equation of motion during close 
approach to a gravitational force center. These techniques reduce or eliminate 
computational and/or analytical problems in calculating those parts of space object 
trajectories [6] . It was shown in previous investigations [7,8,9] that the use of 
regularizing transformations in the formulation of the trajectory optimization problem 
may reduce the computation time and improve the convergence characteristics in 
comparison with unregularized sets of equations. Based on those results, in this study the 
equations for the optimal trajectory of a space vehicle with a continuous low thrust 
propulsion system are derived using regularized variables. The set of variables chosen was 
first described by Sperling [10] and Burdet [11]. The regularization for the state as well 
as for the Lagrangian multiplier equations is obtained by using only the classical 
Sundman [12] time transformation. 

To investigate the numerical behavior of the derived system, example classes are 
chosen: first, numerical calculations of Keplerian orbits and, second, two-dimensional 
earth minimum time escape trajectories starting from various orbits and using various 
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vehicle characteristics. The results are compared in the tirst case with some other problem 
formulations and the analytical solution of Keplerian orbits and, in the second case, with 
the classical unregularized set of Newtonian equations ot motion [6,8,13,14]. The 
comparison indicates that in the cases investigated, a reduction in computation time was 
achieved and that the regularized set of equations of motion is less sensitive to errors in 
the guesses of the unknown boundary conditions. 

A complete analytical solution of the coast arc problem is presented. Some 
remarks on using the new, substituted, redundant variables h and e as independent state 
variables are given. 


PROBLEM FORMULATION 


The equation of motion of a spacecraft with a continuous thrusting propulsion 
system in an inverse-square gravitational force field can be written as 


r 



( 1 ) 


where r is the radius vector from the force center to the vehicle, which is assumed to be 
a point mass; r is the absolute value of the radius vector; p — 7M, where T is the 
universal gravitational constant and M is the mass of the central body; and P is the 
vector of the thrusting acceleration defined as 

P = - u (2) 

m 


where F is the absolute value of the thrust, which is assumed to be constant in this 
study; m = m 0 - j3t, the vehicle mass at time t with 0 as the constant mass flow rate; 
and ii is the unit vector of the unspecified thrust direction. 

The trajectory has to connect the given initial state 


r 0 = f(t 0 ) r 0 = r(t 0 ) 


m 0 = m(t 0 ) 


(3) 


and the final state 
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ff - r(tf) 


ff = r(t f ) 


(4) 


The final time is either predetermined or free. 

Then, the optimal transfer problem considered can be formulated as follows: 
Determine the unspecified thrust history u to minimize any performance ipdex J, a 
function of the boundary values of the state x , where xT = (r, f, m),and the time t. 


J J(x 0 , to, X|-, t^) 


while satisfying the differential equations 


x = f(x, u, t) 


( 6 ) 


for t 0 < t < t f and the boundary conditions given in equations (3) and (4). 

The necessary conditions for an optimal trajectory can be formulated as follows 
(e.g., see Reference 15). After introducing Lagrangian multipliers 


Xj = Aj (t) 


H = constant 


(7) 


the Hamiltonian of the system follows as 


H (x, X, u, ,u) = Yj \ x i * 


( 8 ) 


where f is the constraining equation of the control 


f = (u u) - 1 = 0 


(9) 
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Thus, the Euler-Lagrangian equations are 


— 3H 

3x[ ’ Xi “ " 



( 10 ) 


and the Weierstrass condition is 


H (x 0 , X, u, 0) < H (x 0) X, u 0 , 0) 


(ID 


for each u(t) which fulfills the constraint equation (9) and where x 0 and u 0 describe 
the optimum trajectory. 

Combining equations (1), (2), (8), (10), and (11), the following differential 
equation system is derivable describing an optimal trajectory arc: 


r + 




m IX I 


X + 



3p 


(X • r) _ 
— — r 


( 12 ) 


For a specific mission (e.g., earth escape, orbital transfer, etc.) the additional boundary 
conditions have to be considered. 

Obviously, for impact trajectories r -»• 0, both equations (12) for the state and 
the costate description contain a singularity. Even for close approach to a gravitational 
force center this singularity causes difficulties in numerical calculation of those 
trajectories. From this point of view it might be desirable to find a transformation of 
equations (12) that removes this singularity. In celestial mechanics, regularizing 
transformations are known which remove this singularity within the state equation. 
Applying those regularizing transformations to the optimal trajectory problem shows that 
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some of them will regularize the costate equation, too. In this investigation the 
well-known Sperling-Burdet variables and the classical Sundman time transformation are 
used. In future equations, system (12) is called “Newtonian Formulation of the problem 
in contrast to the regularized form. 


REGULARIZING TRANSFORMATION OF THE 
OPTIMIZATION PROBLEM 


Derivation of the Regularized Equations 

The derivation of regularized differential equations for the state variables and the 
Lagrangian multipliers can be approached either by regularizing the equations of motion 
and then setting down the necessary criteria for optimization in terms ot the new 
independent variable or by setting up the complete optimization problem with ordinary 
time as the independent variable and then applying the regularizing transformation to the 
whole system of resulting equations. In the following derivation, the first approach is 
used. The second approach leads to the same set of equations, as is shown in Appendix 
A, for the second-order system of differential equations. The same connections are 
derived for a differential equation system of the first order by Czuchry and Pitkin [16]. 

Before applying the new independent variables, defined by the classical 
regularizing transformation of Sundman [12] 


ds = r* 1 dt 


(13) 


to the equations of motion, one may write equation (1) as 


r 2 r = - r (r • r) + 2h r - p e + r 2 P 


h = (F-P) (,4) 

I = - [r x (r x P) + P x (r x r)l 
M 


where the following substitution equations for energy h and the Laplacian 
constant e and their time derivatives were used: 
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(15) 


h 


(r • r) m 
2 ‘ r 


e = - - + - [r x (r x r)] 
r n 

Because of equation (13), the transformation rule for a first-order derivative is 


(16) 


d(0 = 1 d(0 
dt r ds 

and a second-order derivative is 

d 2 (-) _ 1 
dt 2 r 3 


d 2 () dr df) 
ds 2 ds ds 


[equation (13a)] 


[equation (13b)] 


Therefrom follows equation (17) for the second time derivative of the radius vector: 


r 


2 f 



(17) 


Combining equations ( 13 ), (14), (17), and the differential equation of the mass 
variation m = j3 and transforming all first time derivatives with the aid of equation (13a) 
gives the complete set of regularized differential equations for the state variables. The 
radius vector equation appears in a second-order form. 

Using the substitution r' = v, the equations can be written in a first-order form: 


r' = v 

v' = 2hr - /re + r 2 P 

h' = (v • P) (18) 
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m' = j3r 
t' = r 


[(18) Concluded] 


Setting up the Hamiltonian of system (18) and deriving its partial derivatives with 
respect to all of the state variables will lead to the following set of Lagrangian multiplier 
equations, the necessary conditions of optimality: 


----- r r 

V - -x v 2h - (A v P)2f - [V X (X e X P) + P X (X e x v)l hi - X m j5 - - X t - 

V = + [ f x ( *e X P) + x e x (f x P)]/m 

x h ' = (X v • r)2 


Xg XyP 


X 


m 



iHpl 


(19) 


X t ' = o 


|H P 1 is the absolute value of the partial derivative of the Hamiltonian H with respect to 
the thrusting acceleration P: 

H p = t 2 X v + A h v + - [r x (v x X e ) + X e ) x (v x ?)] (20) 

M 


The optimal thrust vector program is derived from the condition that the partial 
derivative of the Hamiltonian with respect to the control u, which obeys the constraint 

equation 



(u ■ u) - 1 = 0 


[equation (9)] 


has to be zero. It follows that 


IHpl 


( 21 ) 


The sign in equation (21) follows from examining the second variation of the 
Hamiltonian with respect to u for minimization. 

Equations (2), (18), (19), (20), and (21) give the complete set of equations for an 
optimal trajectory of a thrusting space vehicle in an inverse square force field. The system 
is completely regular for the state description. The equation for A r contains the unit 

vector of the radius r/r which has a well-known limit value for r -*■ 0 (e.g., see 
Reference 17). 

Unfortunately the existence of this limit value, which makes the adjoint system 
regular as well as the state equations, does not prevent numerical difficulties in 
calculating trajectories which contain close approaches to a gravitational force center. For 
calculating plane trajectories the use of polar coordinates is advantageous because the 
adjoint system is completely regular (see following section). The increase of the number 
of differential equations is not an advantage with respect to the integration time. 
Compared to the classical formulation of the trajectory optimization problem [equations 
(12)] in the three-dimensional case, the number of equations increases from 7 to 12 for 
the state description and from 7 to 1 1 for the adjoint system, including the differential 
variation for the mass variation. Those additional equations are one vector equation for 
the Laplacian constant and two scalar equations for the energy and the time, which is 
now a new dependent variable. The adjoint system increases one equation less than the 
state system because the multiplier related to the time t is a constant, as is pointed out 
in Appendix A. If the time would appear explicitly in any of the equations (18) and 
(19), it always could be eliminated by using equation (7). 

To calculate an optimal trajectory, it is only necessary to integrate 22 differential 
equations because of the linear connection between time and the mass variation. The 
equivalent numbers for the problem expressed in the K-S transformation [9] are 11 for 
the state and 10 for the multipliers. However, as shown by Tapley [18], the number of 
differential equations necessary to integrate when solving the optimization problem 
formulated in the K-S transformation is only 18. 

With respect to the number of differential equations for two-dimensional 
trajectory calculations, the use of polar coordinates again may be advantageous because 
the number decreases to seven for the state description and to five for the adjoint system 
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(see following section). The integration interval changes from tQ < t < tf to Sq< s < Sf, 
with Sf unknown even in the case that tf is given because Sf depends on the unknown 
history of the radius of the optimal trajectory [16]. 

In contrast to the classical set of equations for the optimal trajectory, the 
equation of the Lagrangian multiplier associated with the differential equation of the 
vehicle mass X m is not independent. 

In Appendix B, a second-order form of the differential equation tor X v is 
presented. 


Use of Polar Coordinates 

As pointed out in the previous section, the use of polar coordinates especially tor 
calculating two-dimensional trajectories using the regularizing transformation presented 
may be advantageous. 

Transforming equations (18) and (19) in polar coordinates gives the following set 
of equations for the state description: 


v r ' = 2hr + n + r 2 P r 
v ' = r 2 P 

V 

= ^6 

= y l ( r cos 6 ) (22) 

O' = Vg/r 

h' = v r p r + + v 0 P 0 

m' = i3r 
t' = r 


and the following set for the adjoint system: 
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V = • A v 2(h + rP r) ■ 2rP v? - W + VV (r2 COS 0) ' + X 0 v 6 r ~ 2 - \' d 2rP 0 

V = 0 

V = - X^ sin e (r cos 2 #)"’ 

■X Vr X r - Xh P r 

x v ^ = - x h - X^ (r cos 0)-' (23) 

\q = - x h p 0 • V 1 

V = - \ 2r 
K “ ■ ^ 'Hpl 
x t ' = 0 


Obviously this system [equations (22) and (23)1 is in a three-dimensional case not 
regularized. But in two-dimensions, it reduces to 


* = V r 

v r ' = 2hr + ix + r 2 P r 

V = (24) 

h ' = v r P r + W 
m' = (Sr 
t' = r 


and 
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V = - X Vf 2 (h + rP r ) - X Vv? 2r] V - X m0 

V = 0 

X v ;=- X r- X h p r 
= " X h p v? 

K = - ~2 lS P' 

ill m i 

= 0 . 


( 25 ) 


(Hp is defined in previous section.) 

To solve the two-dimensional optimal trajectory problem, the equation system 
necessary to integrate is completely regular. It is only 'necessary to integrate 10 equations 
because the solution is independent of the central angle the connection between mass 
and time is linear, and the Laplacian constant does not appear in the equations o 
motion. To determine the central angle, one additional integration is necessary after the 
optimization problem is solved. 


THE COAST ARC PROBLEM 


Analytical Solution of the Coast Arc Problem 

A coasting arc is defined by setting the disturbing acceleration to zero. From 
equations (18) and (19) the following equations are obtained with P = 0 and m - 0: 

r" - 2h 0 r + (ie 0 = 0 

h' = 0 

e' = 0 (26) 

m' = 0 
f = r 
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and 


X v " - 2h 0 X v = 0 
X h ' = - 2 (X v • f) 

Xg — ^ Xy 

X m ' - 0 
x t ' = 0 . 


( 27 ) 


ho m o are defined by the initial conditions and remain constant over the whole arc. 
System (26), (27) assumes that the coast arc solution is independent of the time t and 
therefrom follows X^ = 0. 

Equations (26) are independent of equations (27). Both sets of equations (26) 
and (27) have an analytical solution. Equations (26) describe Keplerian orbits; they are 
used for numerical investigations in the following section. 

In the case of h 0 < 0. elliptical motion, the solutions of equations (26) and (27) 
are 


r - a cos cos + b sin cos - — ju e 0 

co 2 

t = C 0 + Q s + C 2 cos cos + C 3 sin cos 


(28) 


and 


X v = d cos cos + e sin cos 

r p _ 

X c — — (- e cos cos + d sin cos) + f 

X h = f 0 + f,s + f 2 s 2 + f 3 sin cos + f 4 cos cos + f 5 ssincos + f 6 seoscos 


(29) 


with co 2 - - 2h 0 and a, b, c^, d, e, f, fj constants determinable from the initial 
conditions. 
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The solution of the hyperbolic case h 0 > 0 is similar using hyperbolic instead of 
trigonometric functions. 

For h 0 = 0 parabolic case, the solutions of equations (26) and (27) are 


r = - fj. e 0 s 2 + g 0 s + gj 

(30) 

r s f 

t = J y/ocQ + ttjS + a 2 s 2 + a 3 s 3 + a' 4 s 4 ds 

s o 

and 


\ v = k 0 s + k, 

V = mQ k 0 s 2 + kjS + k 2 ^ 

x h = e 0 + £ t s + e 2 s 2 + e 3 s 3 + V 4 + £ 5 s 5 + e 6 s 6 


(31) 


The constants gg, gj, aj, kg, kj, k 2 , Ej are functions of the initial conditions. 

The solution for t is easily derivable as a sum of elliptical integrals of the first 
and the third kind (e.g., see Reference 19). 

Formulations for the state variable equations of system (26) which allow 
simultaneous solutions of the system for all kinds of orbits may be lound in the 
literature [20] . These formulations are adaptable for the multiplier equations (27), too 
(e.g., see Reference 9). 


Numerical Calculations of Keplerian Orbits 

The numerical calculation of Keplerian orbits is qualified as a test problem for 
the investigation of the numerical stability of the system [14] because a complete 
analytical solution is known. The accuracy of the numerical integration for short and 
long time orbits was compared with those of other regularized sets and the Newtonian 
unregularized set of trajectory equations presented in a paper of Baumgarte [13]. The 
sets compared are 
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1. Set N: The Newtonian equations 


r + 



0 


[equation (1)] 


2. Set ST: A stabilized set (Baumgarte) 


r - - (r • ?')?' + /— in 2 + h ] r = o 

r 2 \2r 2 ) 

t'" + 2ht' = M 


(32) 


3. Set OCS: The oscillator equations (Stiefel, Scheifele K-S Transformation) 



t' = u, 2 + u 2 2 


(33) 


The transformation equations for the state x,, x 2 (the coordinates of the radius vector 
in two dimensions) are x, = u, 2 - u 2 2 , x 2 = 2u!u 2 . 

4. Set REG: Equations (12) with P = 0 


r" - 2hr = - ne 
t' = r 


[equations (26)] 


The unit system for numerical calculations is defined by a = 1 , where a is the semimajor 
axis, and n = 1. The error in the distance due to numerical integration shown in Table 1 
is defined as 


A = [(Ax,) 2 + (Ax 2 ) 2 ] 


1/2 


^ X i x i exact ' x i calculated 


( 34 ) 



TABLE 1. COMPARISON OF THE NUMERICAL INTEGRATION ERRORS 
IN CALCULATION OF KEPLERIAN ORBITS WITH DIFFERENT 
SETS OF EQUATIONS OF MOTION 



2 Revolutions, Starting Point: 
Pericenter 

1 Revolution, Starting 1 
Apocenter 

*oint: 


50 Rev 

olutior 

is 

Eccentricity 

E 

N* 

ST 

osc 

OSC b 

REG 

REG b 

ST 

OSC 

OSC b 

REG 

REG b 

ST 

osc 

REG 

REG b 

0 

7.7' 6 

1.6' 6 

1.7 * 7 

4.8‘ 8 

1.6' 6 

5.8" 8 










0.2 

2.2' 5 

1.2' 5 

2.r 7 

4.8' 8 

US' 6 

S.7' 8 










0.6 

1.2’ 2 

3.3' 6 

3.4- 7 

3.8~ 8 

1.3- 6 

5.0 8 










0.8 

1.6 

4.9' 6 

5.1- 7 

2.8" 8 

9.8" 7 

4.3' 8 










0.9 







2.0' 5 

7.5' 7 

4,7 8 

3.6" 7 

1.9 8 

5.3 4 

8.2" 5 

1 .7’ 5 

8.r 7 

0.95 







l.S' 5 

l.f 6 

4.8 8 

2.6’ 7 

1.8 8 

7.6 4 

1.2 4 

1.2" 5 

8.r 7 

0.99 







l.l 4 

2.6" 6 

4.9 8 

1.2’ 7 

1.6- 8 

5.3' 3 

2.7 4 

5.9’ 6 

8.1 7 

0.999 







7.6' 1 

8.2' 6 

4,9 8 

5.6’ 8 

1.6 8 





. 0.999 999 






- 


5.3' 5 

4.9' 8 

4.3' 8 

1.6" 8 






a. ( )' n means < ) • 10 n 

b. Integration performed with a Runge-Kutta procedure of the fifth-order with step size control 


For the numerical calculations, Runge-Kutta-Fehlberg [3] procedures of the fourth order 
with constant step size (as in Reference 16) and fifth order with step size control were 
used. The high stability of set REG is evident. The constant term pe in equation (26) 
produces stabilization in the long time calculations, even in comparison to the 
similar-looking oscillator-type equations (33). 


MINIMUM TIME EARTH ESCAPE 


Boundary Conditions 

For the minimum time earth escape of the continuous thrusting vehicle 
with j3 the constant mass flow rate, the performance index was formulated as 


J = (-m f ) 


(35) 
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The only boundary condition at the final time is 


h = 0 


( 36 ) 


Therefrom, an auxiliary function 0 is constructed 


4 > = - nrif + 


(37) 


The final values of the Lagrangian multipliers at the final time tf with the 
definition Xj(tf) = Xjf are 


X if = 0 

x mf = -1 (38) 

X hf = v - 

The first equation of set (38) indicates that all final X values, except X m f and X^f are 

zero, v may be expressed in terms of the unknown state variables of the final 
time tf using the condition that the Hamiltonian has to be equal to zero even at the 

final time. Using this condition and applying equations (38) gives 


IjSlrf mf 

v - - 

Flv f l 


(39) 


In deriving equations (38), it was assumed that energy h and Laplacian 
constant e can be considered as independent state variables. Independent is used in the 
sense that h and e are not functions of other state variables, i.e., the partial derivatives 
of <j> [equation (37)] with respect to the radius and velocity are zero. 

The question of the effect on the optimization problem of using additional 
(redundant) variables like h and e are not fully explored here. Some remarks concerning 
this problem are given in Appendix C. 
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Discussion of Numerical Examples 


Equations (2), (18), (19), (20), and (21) and the boundary conditions (3), (36), 
and (38) form the complete set which is to be solved for the two-dimensional optimum 
escape trajectory calculation. The reason that the optimal trajectory is not determined by 
one integration of this set is that the state variables and multipliers are not given as a 
complete set at one definite time t. So the solution is an iteration process starting from 
any initial guess of the missing boundary values at the initial or final time. Usually, an 
initial guess of the Lagrangian multipliers is made. 

In this study approximations of the initial values of the multipliers were 
determined from a backward integration of an escape trajectory using a tangential thrust 
program. This method made use of the knowledge that, for escape, a tangential thrust 
program is close to the optimal steering program. For the integrations, standard 
Runge-Kutta-Fehlberg [3] procedures of different orders as indicated in Tables 1 and 2 
were used. 


TABLE 2. NUMERICAL INTEGRATION CHARACTERISTICS FOR 
DIFFERENT TYPES OF OPTIMAL EARTH ESCAPE TRAJECTORIES 


Trajectory 

Type 

Set of 
Equations 

Coordinate 

System 

Max Error 
Per Step 

Integr. Time 
Per Step (%) 

No. of 
Steps 

Error in 

Quadratic 

Error 0 


‘f 

Spiral 

Newtonian 

Rectangular 

10' 5 

too 

52 

< 10' 7 

< 10' 8 

< 10‘ 12 

'""'3.7 Rev. 

Regularized 

Rectangular 

10' 5 

~ 2130 

37 

< 10" 8 

< 10' 7 

< 10' 13 

t f = 15.45 hr 

Regularized 

Polar 

10 4 

o 

o 

24 

< LO' 7 

< lcf 9 

A 

o 

£ 

Spiral 









~ 240 Rev. 

Regularized 

Polar 

10 4 

~100 

989 

< ID' 6 

< 10' 7 

< 10' 12 

t f = 35 days 









Spiral 









~700 Rev. 

Regularized 

Polar 

10 4 

O 

o 



< 10' 7 

< 10' 12 

165 days 






1 



Highly Ecc. 

Newtonian 

Rectangular 

It)' 5 

100 

26 


< 10" 7 

<10' 13 

t f = 111 hr 

Regularized 

Rectangular 

10' 5 

~20O 

14 


< 10' 6 

<10' 12 

Highly Ecc a 

Regularized 

Rectangular 

10 4 

~100 

132 

O 

V 

A 

O 

CO 

<10' 13 

tf = 56 days 










a. Integration with Runge Kutta 4(5), all other cases with Runge Kutla 7(8) on MSFC’s SDS 930 computer. 

b. Hamiltonian of the system considered. 

c. Sum of the quadratic errors in the boundary conditions at the time t f . 
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The first example calculated is a 15.5 hour earth escape spiral from a low circular 
earth orbit with r = 1.05 times earth radius and an initial thrust to mass ratio of 0.1 
m/sec 2 . This example was also investigated in Reference 14. Although this t/m 0 is 
relatively large, it represents a compromise between a computationally expensive realistic 
trajectory and an inexpensive unrealistic trajectory. The optimal trajectory is shown in 
Figure 1 with the optimal thrusting angle measured against the path tangent. 

Figure 2 shows the relationship between the new independent variable s, the 
regularized time, and the real time t. 




X (1.06 x EARTH RADIUS) 


REGULARIZED TIME 


Figure 1. Optimal low thrust earth 


Figure 2. Real time compared to 


escape spiral. 


regularized time. 


The second example investigated is shown in Figure 3. It is an optimal minimum 
time earth escape starting from a highly eccentric earth orbit near the earth with r 0 = 
1.34 = earth radius. The gravitational force varies until the final time is reached by 
approximately a factor of I0 3 in magnitude; thus, it is expected that for certain error 
bounds the numerical integration step size will change within a wide range. The initial 
thrust to mass ratio of 3.1 O' 3 m/sec 2 is somewhat more realistic than is the spiral escape. 
Figure 3 shows the optimal steering program. Figure 4 shows the relationship between 
the regularized time and real time. The spacecraft in both cases is considered as a point 
mass with a continuous thrusting device. The mass flow rate is constant. All examples are 
calculated in two dimensions. 
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V (1.34 x EARTH RADIUS) 



REGULARIZED TIME 

Table 2 gives some insight into the Figure 4. Real time compared 

integration characteristics of the to regularized time, 

regularized set derived in comparison to 
the Newtonian set of equations ot motion, 

which was also investigated in a first-order form. Some results, not shown in lable 2, 
indicate that the use of a second-order formulation of the Newtonian set of equations 
may save some integration time, it did not improve the convergence behavior of the 
system in those cases investigated. 

The examples for the Newtonian set of equations were calculated in a rectangular 
coordinate system. For the regularized set of equations, rectangular and polar coordinates 
were used. 

The integration time per step is given only in normalized form lor comparison 
with the Newtonian case, which is chosen as 100 percent. The actual values for the SDS 
930 used for the calculations are not representative of modern, faster computers (i.e., 
Univac 1108). The values given are based on a rough calculation from the duration of 
one complete trajectory integration divided by the number of steps. The number ol step 
size changes is not considered but usually there was one additional step each time. 

The doubling in the integration time per step for the regularized set in rectangular 
coordinates in comparison to the Newtonian set is essential due to the increased number 
of differential equations. That also explains the approximate equality ot the integration 
time in the polar coordinates where the number of equations is equal to that of the 
Newtonian formulation. 
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The number of steps for the same error bound required in the boundary 
conditions at the final time is decreased by approximately 30 percent in the 15 hour 
spiral case and in the 111 hour escape case more than 40 percent as compared with the 
Newtonian set and the regularized set calculated in rectangular coordinates. The use of 
polar coordinates for this particular example seems advantageous, even when the number 
of steps is considered. 

Since a closed form solution to the problem considered here does not exist, the 
error generated by the numerical integration process is unknown. The converged 
trajectory is always chosen as the optimal trajectory. 

The error in the Hamiltonian given in Table 2 indicates that all calculations are 
made .within the same error limits and gives an impression of the numerical accuracy 
achieved . 

A certain control of the results is possible by comparing the final time with time 
necessary to achieve zero energy with a tangential thrust steering program. For each 
problem considered, the isolated trajectories calculated with different sets of equations 
and coordinate systems coincide at least within the number of digits indicated by the 
maximum allowable error per integration step. 

The quadratic error given in the last column of Table 2 is defined as the sum of 
the quadratic error of the boundary conditions at the final time. The routine requires 
only that the value be smaller than lCT 12 . The extent to which this criterion is exceeded 
is not controlled. 

To demonstrate the numerical stability of the derived regularized set, three long 
duration optimal escape solutions are shown in Table 2. The 35-day spiral escape 
trajectory starts from a 1.05 earth radius circular orbit, with an initial thrust to mass 
ratio of 10' 3 m/sec 2 . The number of revolutions is more than 235. The 165 day spiral 
escape starts from the same circular earth orbit with an initial thrust to mass ratio of 
510 4 m/sec 2 . The number of revolutions is more than 700. The 56-day escape 
trajectory starts from the highly eccentric orbit. 

Table 3 shows the number of iterations and the number of trajectories calculated 
to isolate the optimal trajectory for all those cases listed in Table 2. The initial guesses of 
the unknown Lagrangian multipliers were in all cases systematically derived from a 
backward calculation starting from the final state, which was obtained from a tangential 
thrust steering program. The deviations of those X values from the finally calculated 
optimal values differ depending on the individual sensitivity of those multipliers with 
respect to a nonoptimal steering program. The deviations range from 10 percent to about 
an order of magnitude and, most importantly, include changes in signs. Table 3 also 
shows the improved behavior of the regularized set in performing those calculations. The 
figures given are possibly not the best values attainable. The influence of different scaling 
systems or optimal selection of sensitivity factors for the optimization technique used is 
not considered for this comparison. 
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table 3. convergence characteristics of different kinds of 

EARTH ESCAPE TRAJECTORIES (Xj CALCULATED BY A 
TANGENTIAL THRUST BACKWARD INTEGRATION) 


Type of 
Escape 

Set of 
Equations 

Coordinate 

System 

Max Error 
Per Step 

Quadratic 

Error 

No. of 
Iterations 

No. of Trajectories 
Calculated 

Spiral 

Newtonian 

Rectangular 


< 10‘ 12 

19 

26 

IS. 45 hr 

Regularized 

Rectangular 


< 10' 12 

6 

16 


Regularized 

Polar 


< 10' 13 

4 

11 

Spiral 
35 days 

Regularized 

Polar 

10 4 

< 10' 12 

30 

37 

Spiral 
165 days 

Regularized 

Polar 

10 4 

< to' 12 

18 

25 

Highly Ecc. 

Newtonian 

Rectangular 

10' 5 

n 

o 

V 

29 

36 

Ill hr 

Regularized 

Rectangular 

10 4 

< 10- 15 

6 

16 

Highly Ecc. 
56 days 

Regularized 

Rectangular 

10 4 

< 10' 13 

13 

23 


The investigator was able to calculate the long duration cases only with the aid 
of the regularized set of equations of motion. 

To investigate the sensitivity of the different sets of equations with respect to 
deviations of the initial X values from the optimal values in Tables 4 and 5, the results 
of the non-realistic systematic changes in the initial X’s are given for the shorter escape 
cases. It is found that for small deviations due to the smaller number of differential 
equations, the Newtonian set may require less trajectory calculations than the regularized 
set; but for greater deviations the regularized set is obviously less sensitive. The remark 
“non-converged” indicates that after about 50 iterations, no improvement of the isolation 
procedure was noticed. Also, for the figures in Tables 4 and 5, no effort was made to 
derive the smallest possible number of iterations. For the sensitivity analysis, the error in 
the final time was always considered to be zero. 
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TABLE 4. INITIAL ERROR INFLUENCE ON THE CONVERGENCE 
CHARACTERISTICS OF THE 15.45 HR EARTH ESCAPE SPIRAL 


\ / initial X \ 

Newtonian Set of Equations 
in Rectangular Coordinates 

Regularized Set in 
Rectangular Coordinates 

\>pt ^optimal X J 

N a 

n c 

N 

n 

0.01 

0.1 

Not Conv.k 
Not Conv. 


Not Conv. 
21 

31 

0.5 

17 

24 

8 . 

18 

0.8 

7 

14 

5 

15 

1.0 

0 

1 

0 

1 

1.2 

7 

14 

6 

16 

2.0 

16 

23 

16 

26 

10.0 

19 

26 

20 

30 

100.0 

Not Conv. 


30 

40 


a. N - Number of iterations. 

b. Not converged within 50 iterations. 

c. n - Number of trajectories calculated to Find the optimal solution. 


TABLE 5. INITIAL ERROR INFLUENCE ON THE CONVERGENCE 
CHARACTERISTICS OF THE 1 1 1 HR EARTH ESCAPE TRAJECTORY 


^i / initial X \ 

Newtonian Set of Equations 
in Rectangular Coordinates 

Regularized Set in 
Rectangular Coordinates 

\>pt \ optimal \J 

N a 


N 

n 

0.1 

Not Conv. 


20 

1 

0.5 

35 

42 

10 

1 

0.8 

38 

45 

C 


1.0 

0 

1 

0 

i 

1.2 

20 

27 



1.5 

38 

45 

8 

18 

10.0 

69 

76 

20 

30 


a. N - Number of iterations. 

b. n - Number of trajectories calculated to find the optimal solution. 

c. Not calculated. 
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CONCLUSIONS 


Based on the results shown in this paper, the following conclusions may be 

drawn. 


For the examples calculated, the integration time per trajectory necessary in the 
Newtonian and regularized formulations of the equations of motion are approximately 
equal. The gain due to the smaller number of steps using the regularized set is lost 
because of higher integration time per step which results from the use of a greater 
number of differential equations. 

In the case of good initial guesses of the missing boundary values (i.e., the 
Lagrangian multipliers at the time t = 0), the smaller number of differential equations in 
the Newtonian formulation is even advantageous with respect to the numerical 
optimization method used. Only slight differences between both formulations were 
found. 


On the other hand, a remarkable gain in computation time is achievable for the 
trajectory types investigated by using the regularized set in the more realistic cases of bad 
initial guesses for the Lagrangian multipliers. The smaller sensitivity of the regularized set 
leads more rapidly to the converged trajectory. In the case of extremely bad initial 
values, the regularized set gives the hope that convergence may be achieved where it 
otherwise may have been impossible. 

In the long duration cases presented, convergence was achieved without use of 
aids other than those mentioned above and without major difficulties. 

George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 

Marshall Space Flight Center, Alabama, October 1973 
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APPENDIX A 


AN ALTERNATE WAY TO DERIVE 
THE REGULARIZED EQUATIONS 


Using the equation of motion [equation (1)] with the substitutions for 
energy h and Laplacian constant e that were made in equations (15) and (16) and the 
additional substitution of the velocity v = r, one finds the first-order system of the 
equations of motion: 


r 


v 


1 

v = - - [v (r -v) - 2h r + ge - r 2 P] 
r* 

li = (v - P) (A-l) 

A 1 - 
e - - (v x (r x P) + P x (r x v)) 

P 

m = 0 


If one starts setting up the optimization problem with the system (A-l) and using 
the time transformation dt = rds, one will not tind the same set of Lagrangian multiplier 
equations (19). In particular, the system derived in such a manner will be more complex 
and not regularized, and the Lagrangian multiplier equation related to the differential 
equation of the vehicle mass will stay independent as in the classical case. However, 
starting from the second-order equation for the radius vector 


r 2 r = - r (r ■ ?) + 2h r - pe + r 2 P 


(A-2) 


and considering the tact that the regularized velocity using the time transformation dt = 
rds will be 


(A-3) 


24 



the following first-order system will be derived : 


r = v/r 

v = (2h r - jue + r 2 P)/r 

li = (v ■ P)/r (A-4) 

~ [v x (r x P) + P x (r x v)] 

e = - 

r 

m = j3 


Starting with set (A-4) to set up the optimization problem gives the same Lagrangian 
multiplier equations as in the first approach described above. 

The relation between S H, the Hamiltonian of the regularized set of state 
equations (12), and t H, the Hamiltonian of the unregularized set (A-4), is 


s H = r l H 


(A-5) 


Building the partial derivatives gives the connection between the regularized and 
unregularized sets of multiplier equations: 


V 



% 


- (r l H) = 
9r 


a s H 

"a r 


- \ r 
X h = X h r 
x e ” X e r 
x m _ X m r 


(A-6) 
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The second part of the first equation in set (A-6) appears because of the defining 
equation of the new independent variable. The additional multiplier X t follows for the 

autonomous system (18) and, as pointed out above, considering the necessary free 
boundary condition Sf as the constant, 

X t = - *H . (A-7) 

is the value of the Hamiltonian of the unregularized system (A-4) (see Reference 16). 
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APPENDIX B 


SECOND-ORDER FORM OF LAGRANGIAN 
MULTIPLIER EQUATION X v 


The Lagrangian multiplier equations associated with the radius vector* and 
velocity, equations (19), may be written as 


X 


r 


- 2hX v 


3Hp _ f r 

"aT p + v> ; + ; 


(B-l) 


and 


- - V 



(B-2) 


Finally, the second-order multiplier equation is 



9Hp _ 

P + 

3r 



- + x t> 


r 

r 


(B-3) 


and the second-order radius equation remains 

r" - 2hf = - + r 2 P . (B ' 4) 

The other equations for h', t, t', m\X h ', X e \ X t ', X m ' are unchanged [see equations (18) 
and (19)] . 

In comparison to the “Newtonian” Formulation, equations (12), the second-order 
multiplier equation for X v contains first-order derivatives on the right-hand side. 
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APPENDIX C 


INDEPENDENCE OF 

NEW SUBSTITUTED VARIABLES LIKE h AND e 


Instead ot systems (12) and (18), consider a simple optimization problem in 
general formulation. 

Given is a differential equation system 

* = f(x, u, t) (C-l) 

and a definition equation 


y = y(x, u, t) 


(C-2) 


Substituting y [equation (C-2)] into equation (C-l), the system of differential equations 
can be written as 


x = f(x, y, u, t) 
y = g(x,y, t) 


(C-3) 


By intioducing Lagrangian multipliers x X, ^X, ^X, where the left-hand superscript 
denotes to which of the state variables the adjoint variable is related, a function H (the 
Hamiltonian) can be defined for system (C-l): 


X H = x X T f 


(C-4) 


and for system (C-3): 


xy H = *x T f + yx T g 


(C-5) 


28 



(For distinction in setting down X ^H for system (C-3) k is used as a superscript 
for X related to x instead of x as used for equation (34). 

Using equation (C~4), the following differential equations for the 
multipliers X A may be derived: 

x X = - X H X = x X T (f x + f y y x ) ( c ‘ 6 > 


where f is used as defined within the first equation of system (C-3) and y defined by 
equation (C-2). The right-hand side subscripts denote partial derivatives. 

From equation (C-5), the following differential equations are derived. 


= - xy H x = - “X T f x - y X T g x 
y X = - xy H y = - K X T f y • 


(C-7) 


In order to show that system (C-6) is equivalent to system (C-7) while under 
consideration of 


g x = (y) x = (y x ) > 


(C-8) 


the first equation of system (C-7) is 

K X = - K X T f x - y X(y x ) + y X T y x - y X T y x 


(C-9) 


and the second equation of system (C-7) is 


K X = - K X T (f x + f y y x ) - ( y X T y x ) 


(C-10) 



Integrating the system (C-10) over t shows that the second part of the right-hand side of 
equation (C-10) is dependent on the boundary conditions only 


t f 

f ( y Xy x ) dt = [ y Xy x ] Q 
0 


(C-ll) 


Thus, it follows that systems (C-7) and (C-6) are equivalent. They differ by an additive 
constant only. 

The derivation of the equivalence of both Lagrangian multiplier systems shows 
that in both cases the same optimal arc will be derived using either this or the other 
formulation. It does show nothing about the influence on the procedure of calculating 
those optimum solutions carrying more, redundant, variables (as in this case y). 
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